#include<iostream>
#include<fstream>
extern "C"{
#include<math.h>
}
#include "IMParton.h"
using namespace std;

IMParton::IMParton(unsigned int Z_temp, unsigned int A_temp):Parton(Z_temp, A_temp)
{
	cout<<"    IMParton version - 0.0"<<endl;
}

IMParton::~IMParton(void)
{
}

double IMParton::getPDF(int Iparton, double x, double Q2) const
{
	if(Iparton<0 || Iparton>5)
	{
		cout<<"****Unknown Iparton value."<<endl;
		cout<<"****Iparton must be one of these: 0, 1, 2, 3, 4, 5."<<endl;
		return 0;
	}
	else
	{
		double lnx, lnQ2;
		int i=(int)(lnx=log(x*1e6)/lnxstep);
		int j=(int)(lnQ2=log(Q2*8)/lnQ2step);
		double g0[3], g1[3], g2[3], g[3]={0};
		if(i<0)i=0;
		if(i>(xMax-3))i=xMax-3;
		if(j<0)j=0;
		if(j>29)j=29;
		//avoid log(1-x) calculation in below algorithm
		if(x>0.999999);
		//if x>0.5, we use A(1-x)^B to do the interpolation
		else if(x>0.5)
		{
			double vln1_x[2]={log(1-exp(log(1e-6)+i*lnxstep)),log(1-exp(log(1e-6)+(i+1)*lnxstep))};
			g0[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g0[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			j++;
			g1[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g1[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			j++;
			g2[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g2[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			g[0]=exp(fitLinear(log(1-x),vln1_x,g0));
			g[1]=exp(fitLinear(log(1-x),vln1_x,g1));
			g[2]=exp(fitLinear(log(1-x),vln1_x,g2));
		}
		//if x<1e-5, we use A*x^B to do the interpolation
		//for valance quark, B>0; for gluon and sea quark, B<0
		else if(x<1e-5)
		{
			double vlnx[2]={i,i+1};
			g0[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g0[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			j++;
			g1[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g1[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			j++;
			g2[0]=log(grid[(xMax*j+i)*6+Iparton]);
			g2[1]=log(grid[(xMax*j+i+1)*6+Iparton]);
			g[0]=exp(fitLinear(lnx,vlnx,g0));
			g[1]=exp(fitLinear(lnx,vlnx,g1));
			g[2]=exp(fitLinear(lnx,vlnx,g2));
		}
		//we use quadratic interpolation method for other situations
		else
		{
			double vlnx[3]={i,i+1,i+2};
			g0[0]=grid[(xMax*j+i)*6+Iparton];
			g0[1]=grid[(xMax*j+i+1)*6+Iparton];
			g0[2]=grid[(xMax*j+i+2)*6+Iparton];
			j++;
			g1[0]=grid[(xMax*j+i)*6+Iparton];
			g1[1]=grid[(xMax*j+i+1)*6+Iparton];
			g1[2]=grid[(xMax*j+i+2)*6+Iparton];
			j++;
			g2[0]=grid[(xMax*j+i)*6+Iparton];
			g2[1]=grid[(xMax*j+i+1)*6+Iparton];
			g2[2]=grid[(xMax*j+i+2)*6+Iparton];
			g[0]=fitQuadratic(lnx,vlnx,g0);
			g[1]=fitQuadratic(lnx,vlnx,g1);
			g[2]=fitQuadratic(lnx,vlnx,g2);
		}
		//if Q2>1, we do the interpolation to the variable ln(Q^2)
		if(Q2>1)
		{
			double vlnQ2[3]={j-2,j-1,j};
			return fitQuadratic(lnQ2,vlnQ2,g);
		}
		//if Q2<1, we do the interpolation to the variable Q^2
		else 
		{
			double vQ2[3]={0.125*pow(2,j-2),0.125*pow(2,j-1),0.125*pow(2,j)};
			return fitQuadratic(Q2,vQ2,g);
		}	
	}
}


